##########################################################################################

library(dplyr)
library(ggplot2)
library(data.table)
library(RColorBrewer)
library(optparse)

##########################################################################################

option_list <- list(
    make_option(c("--maf_im_file"), type = "character") ,
    make_option(c("--maf_im_msi_file"), type = "character") ,
    make_option(c("--images_path"), type = "character") ,
    make_option(c("--info_file"), type = "character")
)

if(1!=1){
    
    work_dir <- "~/20220915_gastric_multiple/dna_combinePublic/"
    maf_im_file <- paste(work_dir,"/maf_public/All_use.IM.maf",sep="")
    maf_im_msi_file <- paste(work_dir,"/maf_public/All_use.IM.msi.maf",sep="")
	images_path <- paste(work_dir,"/images/mutRate",sep="")
	info_file <- paste(work_dir,"/baseTable/STAD_Info.addBurden.MSI_MSS.addCNVType.tsv",sep="")

}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

maf_im_file <- opt$maf_im_file
info_file <- opt$info_file
images_path <- opt$images_path
maf_im_msi_file <- opt$maf_im_msi_file

dir.create(images_path , recursive = T)

###########################################################################################

info <- data.frame(fread(info_file))
dat_maf_im <- data.frame(fread( maf_im_file ))
dat_maf_im_msi <- data.frame(fread( maf_im_msi_file ))

###########################################################################################

Variant_Type <- c("Missense_Mutation","Nonsense_Mutation","Frame_Shift_Ins","Frame_Shift_Del","In_Frame_Ins","In_Frame_Del","Splice_Site","Nonstop_Mutation")

info$Molecular.subtype <- info$TCGA_Class

###########################################################################################
## 肠化样本
maf <- dat_maf_im
maf_use <- data.frame(
    Hugo_Symbol = maf$Hugo_Symbol, 
    Chromosome = maf$Chromosome , Start_Position =  maf$Start_position , 
    Reference_Allele = maf$Reference_Allele , Tumor_Seq_Allele2 = maf$Tumor_Seq_Allele2 , 
    Variant_Classification = maf$Variant_Classification , Tumor = maf$Tumor_Sample_Barcode , From = maf$From )

maf_use <- subset( maf_use , Variant_Classification %in% Variant_Type )
maf_use <- merge( maf_use , info[,c("Patient" , "Class" , "Molecular.subtype")] , by.x = "Tumor" , by.y = "Patient")
maf_im <- maf_use
maf_im$Class <- "IM"

maf <- dat_maf_im_msi
maf_use <- data.frame(
    Hugo_Symbol = maf$Hugo_Symbol, 
    Chromosome = maf$Chromosome , Start_Position =  maf$Start_position , 
    Reference_Allele = maf$Reference_Allele , Tumor_Seq_Allele2 = maf$Tumor_Seq_Allele2 , 
    Variant_Classification = maf$Variant_Classification , Tumor = maf$Tumor_Sample_Barcode , From = maf$From )

maf_use <- subset( maf_use , Variant_Classification %in% Variant_Type )
maf_use <- merge( maf_use , info[,c("Patient" , "Class" , "Molecular.subtype")] , by.x = "Tumor" , by.y = "Patient")
maf_im_msi <- maf_use
maf_im_msi$Class <- "IM"

###########################################################################################
## 总的，不同来源的
## IGC和DGC分开的，合并的
maf_use <- rbind(maf_im , maf_im_msi)

###########################################################################################

msi_sample <- unique(subset( info , TCGA_Class == "MSI" )$Patient)
cin_sample <- unique(subset( info , TCGA_Class == "CIN" )$Patient)
gs_sample <- unique(subset( info , TCGA_Class == "GS" )$Patient)

###########################################################################################

info_im <- subset( info , Class == "IM" & Type != "IM + IGC + DGC" ) 
info_im$Molecular.subtype <- ifelse( info_im$Patient %in% msi_sample , "MSI" , "" )
info_im$Molecular.subtype <- ifelse( info_im$Patient %in% cin_sample , "CIN" , info_im$Molecular.subtype )
info_im$Molecular.subtype <- ifelse( info_im$Patient %in% gs_sample , "GS" , info_im$Molecular.subtype )

info_use <- unique(info_im[,c("Patient" , "Type" , "Molecular.subtype" )])
info_use1 <- info_use
info_use1$Type <- "All"

info_use <- rbind( info_use1 , info_use )

###########################################################################################
## 按照病理类型计算突变率
## 以人为单位
maf_use1 <- merge( maf_use , info_use[,c("Patient" , "Type")] , by.x = "Tumor" , by.y = "Patient" )

mut_rate <- maf_use1 %>%
group_by( Hugo_Symbol , Molecular.subtype , Type ) %>%
summarize( MutNum = length(unique(Tumor )) )
mut_rate$id <- paste0( mut_rate$Type , ":" , mut_rate$Molecular.subtype )

class_num <- info_use %>%
group_by( Type , Molecular.subtype) %>%
summarize( SampleNum = length(unique(Patient)) )
class_num$id <- paste0( class_num$Type , ":" , class_num$Molecular.subtype )

mut_rate <- merge( mut_rate , class_num[,c("id" , "SampleNum")] , by = "id" )
mut_rate$MutRate <- mut_rate$MutNum/mut_rate$SampleNum

out_name <- paste0( images_path , "/MutRate.molType.IM.tsv" )
write.table( mut_rate , out_name , row.names = F , quote = F , sep = "\t" )